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Abstract. Different schemes for the treatment of long-ranged electrostatic interac- 
tions will be examined for water simulations using the polarizable fluctuating charge 
potential. Several different methods are compared, including Ewald sums, potential 
shifting, spherical truncation and reaction field corrections. For liquid water, prop- 
erties such as the energy, pressure, dynamics and structure are more sensitive to the 
treatment of the long-ranged interactions with polarizable than with non-polarizable 
potentials. 



INTRODUCTION 

Simulating systems with long-ranged electrostatic interactions using periodic 
boundary conditions requires a treatment of the interactions beyond the central 
simulation box, either using Ewald sums, truncation or modification of the po- 
tential. The importance of the proper treatment of long-ranged forces has been 
demonstrated for many systems, including pure water [1-5], hydrophobic aggrega- 
tion [6], ionic solutions [3], and proteins and peptides [7-11]. Nevertheless, many, 
if not most, simulations of aqueous systems are done with periodic boundary con- 
ditions but without using Ewald. There are two main motivations for avoiding 
Ewald. First, many seek to avoid the additional computational time that evaluat- 
ing the Ewald sums requires, although algorithms such as the cell multipole [12] and 
particle mesh Ewald [13,14] have made Ewald efficient for large systems. Second, 
many seek to use periodic boundary conditions to avoid edge effects but eliminate 
Ewald in the hope that having no direct interactions between periodic images better 



^) The content of this publication does not necessarily refiect the views or policies 
of the Department of Health and Human Services, nor does mention of trade names, 
commercial products, or organization imply endorsement by the U.S. Government. 



represents infinite dilution for aqueous solutes. As stated by Allen and Tildesley, 
Ewald methods "will tend to overemphasize the periodic nature of the model fiuid 

[15]." An alternative to Ewald and simple truncation are reaction field methods in 
which the volume outside the cut-off distance is treated as a dielectric continuum 
[15,16]. The infiuence of long-ranged interactions may be even more significant for 
polarizable systems since the Coulomb or dipole-dipole interactions will couple to 
the polarizable charges or dipoles. The recent widespread use of polarizable water 
models, sometimes without using Ewald [17-27], suggests that an examination of 
the effects of truncation on polarizable systems is necessary. The polarizable ffuctu- 
ating charge model for water will be used [28] to examine how different truncation 
schemes and Ewald influence the structure and dynamics of pure water. 
Six different simulation methods will be used. 

A. Ewald sums. 

B. Scaling by the complementary error function, Ss(r)=erfc(Ar). This is simply 
the real space part of Ewald and A is set equal to 5/L as is fairly standard [29]. 

C. Scahng by erfc(Ar), with the addition of the Ewald self-term and mean-fleld 
approximation for Fourier space term. 

D. Scahng by 



Loncharich and Brooks used n = 2 in Equation 1 [7]. However, References [3] 
and [4] found that n = 1 works better than n = 2 for pure water. The present 
simulations will use n = 1 and rcnt=L/2, where L is the box length. 

E. Nearest image truncation. Coulombic interactions between two-molecules are 
included only if the distance between center-of- masses is less than a cut-off distance, 
taken to be L/2. 

F. Reaction field correction to truncation. The Coulomb interaction becomes 



where e^^^ is the dielectric constant of the reaction field, Vcut is the cut-off radius and 
rjiF is the radius of the reaction field. Hiinenberger and van Gunsteren have found 
that TRF—Tcut is optimal [5] and so the present simulations will use rRF—f'cut—^/'^ 
and eRF is set equal to 79. With r-RF—rcut and erf much larger than 1, the reaction 
field interaction becomes a scaling function acting on the Coulomb interaction equal 



The different methods are plotted on Figure 1, which compares S(r)/r with the 
bare Coulomb interaction. Notice that the functions SB(r), S£)(r) and the reaction 
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field interaction are similar and modify the potential over its whole range. Notice 
also that at the cut-off, at 10, the Coulomb interaction is 0.1, only 20% lower than 
its value at hydrogen bonding distances of 2 and is still far from zero. Cut-off 
lengths are usually in the range of 8 to 10 A. 

The periodic boundary conditions can either be applied to atoms individually or 
by groups. If applied by groups (for example, by molecules or residues), then the 
entire group is moved together so that the nearest image between group centers 
(for example, the center of mass) is used. If applied by atoms, then each atom 
is moved individually. The second method can create large dipoles by splitting 
molecules across a box length. On the other hand, applying the periodic boundary 
conditions to groups can suddenly shift entire molecules across the box, resulting 
in discontinuities. If the length scale of a molecule or residue (the distance from 
the group center to the outermost atom) is R then the distance between two atoms 
can suddenly change from |r+R| to |r-R|. These problems will not be present if 
the potential equals zero at L/2-R. For the methods in which the potential slowly 
approaches zero at L/2 (methods A-D and F), atom based periodic boundary 
conditions are used, since this was found to conserve energy better. (Although 
the differences between atom and group based periodic boundary conditions are 
not large for the Ewald simulation.) For nearest image truncation with a spherical 
cut-off, the treatment of nearest image makes a large difference and using atom- 
based periodic boundary conditions leads to great instabilities in the system, so a 
molecular based method was used. 

The effects of the various treatments of the long-range forces on the energy, pres- 
sure, dipole moment and pair correlation functions will be examined. In addition, 
the dependence of dynamical properties in terms of the translational and rotational 
diffusion constants will be examined. Convergence of these properties will be exam- 
ined for two different system sizes: 256 water molecules (corresponding to L=19.7A 
at a density of 1 g/cm^) and 512 molecules (L=24.8A). 

POTENTIAL MODEL 

The fluctuating charge (FQ) model is a polarizable potential model in which the 
partial charges on atomic sites are treated as variables which respond to changes 
in their environments [28]. The model gives accurate predictions for liquid state 
properties. At T=298 and P=l atm, the dielectric constant is 79. The TIP4P-FQ 
model uses the geometry of the TIP4P water model and includes Lennard- Jones 
interactions between oxygen sites and three charge sites: two on the hydrogen 
atoms and one on the M-site 0.15A from the oxygen atom [30]. The FQ model has 
additional interactions between charge sites on the same molecule. The energy is 

E = Elj + Ecoulomb + Epol (3) 

where E^j is the sum of all Lennard- Jones interactions between oxygen sites, 
Ecouiomb is the sum of all Coulomb interactions between different molecules 



-'-^Coulomb (4) 

where Qjc is the charge of atom a on molecule i and Epo; is the difference in the 
molecular energy between the liquid and gas-phase, 

Epol = X] [ X] y^aQia + 9 X! X! QiaQipJafiir ia,il}) " Egp j (5) 

where Egp is the gas-phase energy, x|J^ is the MuUiken electronegativity of the iso- 
lated atom and Jai3{fia,ii3) is the intramolecular interaction. The gas-phase energy 
is the minimized energy of the isolated molecule. With Ewald, ^couimnh becomes 

Ecoulomb = ^^Y1 QiaQjP e.^ f c{\r p) / T ia,j (3 (6) 

and there are two additional energy terms: the Fourier space term 

(7) 
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and the self-term, which corrects for including i = j terms in the same box in Ei?5, 

Eself = H QiaQm (^rfiXria,ip)/ria,il3- (8) 

where erf(x) is the error function [15]. The screening parameter A was set equal 
to 5/L, 256 lattice vectors were used in the Fourier space sum and conducting 
boundary conditions were used. There is a self energy for the reaction field method 
also, which is given by [5] 

Eself = 9 QiaQip II ■ 

^ i a /3 ^^RF + -I- Tj^p 

Model C consists of the Coulomb energy from Equation 6 and the self term (Equa- 
tion 8) plus an approximate Fourier space term 

Efs = \y.Y.Q^o.{€') (10) 



with 

G^O ^ i j 13 



which was calculated from the Ewald simulation. These terms are constants for 
rigid molecules but by coupling to the charges they can influence the dynamics. For 
256 water molecules and L=19.7A, {(f)^^) =0.0125 kcal/mol/e and =-0.0280 
kcal/mol/e. 

The set of charges which minimize Equation 3 are the ground state charges, 
subject to a charge neutrality constraint on each molecule. Rather that solving 
for the charges exactly at each time step, the method treats them as dynamical 
variables, which are propagated in an extended Lagrangian formalism at a low 
temperature so as to remain near the potential energy minimum [28]. The extended 
Lagangian method introduces a new complication when using cut-offs. If the cut- 
offs introduce discontinuities in the potentials (as method E does) then propagating 
the charge degree-of- freedom becomes more difficult, because the forces on the 
charges are discontinuous. For the case of method E, in order to prevent the charge 
degrees-of-freedom from getting too hot and drifting away from the potential energy 
minimum, the exact set of charges was found every 1000 times steps. In all other 
cases, the charge temperature remains under 10 Kelvin for a 100 ps simulation. 
In addition, we found that the discontinuities in spherical truncation caused the 
system to gradually heat up (at about 3 K/ps) so these simulations were all done 
with a Nose-Hoover temperature bath with a mass for the Nose variable equal to 2.0 
kcal/mol psec^ [31,32]. All other simulations were done at constant E,V,N, except 
for equilibration (at T,V,N) and as noted. The simulations were done with a 1 fs 
time step and used SHAKE to enforce bond constraints [33]. The Lennard- Jones 
interactions were calculated only between the nearest periodic images. This too 
introduces some discontinuities into the energy and forces and sometimes switching 
functions are used for the Lennard- Jones interactions. However, for a box length 
of 20A, the TIP4P-FQ Lennard- Jones interaction at half the box length is only 
1x10"'^ kcal/mol. The data presented in the next section is from four separate 100 
ps runs. 

RESULTS 

Energy. The energies for the different simulations are hsted in Table 1. The paren- 
theses give 95% confidence intervals. The Ewald results (A) do not show a system 
size dependence for the total energy. The individual contributions (Ecouiomb, Efs 
and ^seif) should be system size dependent, since they use a value of A dependent 
on the box size, but the sum of these three terms should be size independent. The 
largest component of the Ewald terms is Ecouiomb- The other terms Eps and E^eif 
make much smaller contributions. However, removing these terms, which would 
correspond to using a complementary error function shifting function (method B), 
gives a much different energy, which is higher by 1 kcal/mol for the 256 molecule 
system. For the larger system with a longer r^ut the energy is closer to the Ewald 
result. If the self-term and a mean-field estimate of E^^^ (see Equations 10 and 11) 
are added back in (model C), the results are improved considerably. These results 



are essentially indistinguishable from the Ewald results. 

The use of the shifting function (model D) gives energies similar to shifting by 
the complementary error function and the energies show a strong dependence on 
cut-off length. Other studies with non-polarizable water potentials have also found 
that using the shifted potential, S/)(r), leads to an increase in the energy. For 
the SPG potential with rc„t=9.3A, the energy is 0.4 kcal/mol higher [4] and for 
the TIPS potential with rc„j=8.0A, the energy is 0.6 kcal/mol higher [3]. For the 
polarizable model with a similar cut-off distance (for 256 molecules, rc„f=9.85A) 
the difference in the energy is greater. Nearest image truncation (model E), on 
the other hand, overestimates the energy by almost a half a kcal/mol and does not 
improve with system size. For a non-polarizable models, the energy with spherical 
truncation also does not show much of a dependence on cut-off length [2,34] and 
overestimates the energy, but only by 0.1 kcal/mol [5]. Once again the differences 
between Ewald and other treatments is greater for the polarizable model. Another 
study using spherical truncation with rc„i=10.5A and the SPC-FQ model also finds 
a lower energy (-11.5 kcal/mol) [27] than the Ewald result (-9.9 kcal/mol) [28]. For 
the RPOL model of water, which treats polarizability using point inducible dipoles, 
spherical truncation (with rcut=9A) only slightly overestimates the energy [35,36]. 
The cut-offs apparently are more severe for the charge-charge (1/r) interactions 
than for the dipole-dipole (1/r^) interactions. In fact, the cut-offs introduce less 
errors for the RPOL model than for non-polarizable models. 

The reaction field method (model F) gives very good agreement with the Ewald 
results for the energy. This agreement, and the improvement over the shifting 
function (method D), is remarkable considering the similarity of the treatment of 
the Coulomb interaction (see Figure 1). The reaction field interaction (Equation 2) 



TABLE 1. Total potential energy, divided by the number of molecules, and energy 
components for the various treatments of long-ranged electrostatics for two different 
sized systems, in kcal/mol. 





number of 
molecules 


Etot 


Elj 


Ecoulomb 


Efs 


Eself 


Epoi 


A 


256 


-9.86(5) 


2.30(7) 


-17.2(2) 


0.018(1) 


-0.659(4) 


5.6(1) 




512 


-9.85(5) 


2.24(5) 


-17.4(2) 


0.009(1) 


-0.321(1) 


5.6(1) 


B 


256 


-8.85(9) 


1.59(6) 


-14.6(2) 






4.2(1) 




512 


-9.4(1) 


1.90(9) 


-16.1(3) 






4.9(1) 


C 


256 


-9.84(3) 


2.30(4) 


-17.1(1) 


0.025(1) 


-0.658(2) 


5.6(1) 


D 


256 


-9.09(4) 


1.59(3) 


-15.0(1) 






4.4(1) 




512 


-9.34(9) 


1.78(8) 


-15.8(3) 






4.7(1) 


E 


256 


-10.21(3) 


2.5(1) 


-18.9(3) 






6.2(1) 




512 


-10.20(3) 


2.49(2) 


-18.8(1) 






6.1(1) 


F 


256 


-9.92(8) 


2.20(9) 


-17.7(3) 




-0.051(1) 


5.6(1) 




512 


-9.94(5) 


2.31(5) 


-17.9(2) 




-0.026(1) 


5.7(1) 



experiment -9.9"* 
a. Reference [30] 



TABLE 2. Total pressure and pressure components, in kbar. 



number or 


T> 

i^tot 




1 ) 

Coulomb 




molecules 










A 256 


0.02(4) 


50.7(7) 


-52.0(7) 


-0.038(3) 


512 


0.04(5) 


50.2(5) 


-51.6(5) 


-0.015(2) 


B 256 


0.5(2) 


43.3(6) 


-44.2(1) 




CIO 


0.3(1) 


46.7(9) 


-47.8(9) 




L z5d 


0.15(d) 


50.7(4) 


-51.9(4) 




D 256 


-1.5(1) 


43.2(3) 


-46.2(4) 




512 


-1.3(1) 


45.4(9) 


-48.1(9) 




E 256 


-0.7(2) 


53(1) 


-55.4(9) 




512 


-0.77(6) 


52.8(2) 


-55.0(3) 




F 256 


-2.9(1) 


49.6(9) 


-53.9(9) 




512 


-2.4(1) 


51.0(5) 


-54.7(6) 




experiment 


0.001 









is greater than SD(r)/r so the electrostatic interactions are stronger, which leads 
to the lower energy. In addition, there is the self-term (Equation 9) which acts 
to slightly increase the magnitude of the charges and also lowers the energy. For 
non-polarizable water potentials, the energy using the reaction field method also 
agrees well with the energy using Ewald [5] . 

Pressure. The pressure is a balance of repulsive and attractive forces and so is 
sensitive to the treatment of electrostatic interactions (Table 2). Listed are the 
total pressure, plus the contributions to the virial from the different interactions. 
The components do not add up to the total because there is the additional ideal 
gas part. Only the Ewald method gives the correct pressure. The reaction field 
method, which gave an accurate energy, does not do well for the pressure. The 
differences in the pressure between Ewald and other methods are more substantial 
for the polarizable model than for non-polarizable models [4,5]. Like the energy, 
the pressure with Ewald is size independent. A previous simulation using Ewald 
with a non-polarizable potential found that the energy and pressure shows no size 
dependence for systems of 64 or more molecules [37] . 

Dipole moment. The dipole moment, which in all cases is enhanced relative to 
the gas-phase value of 1.85 Dcbye, correlates very well with the energy. Treatments 
which give an accurate energies (models C and F) also have the same dipole moment 
as the Ewald method. Those with a larger dipole moment (model E) have a lower 
energy and those with a smaller dipole moment (models B and D) have a higher 
energy. It is the sensitivity of the charges that makes the proper treatment of 
long-ranged electrostatics more important for polarizable models. The method C, 
which did not work well, is improved considerably just by adding the constant 
terms which couple to the charges to give the right dipole moment. 

Dynamical properties. Also listed in Table 3 are the translational diffusion con- 



stant and the NMR relaxation time, tnmr, which gives the time scale for rotations 
around the axis connecting the hydrogen atoms [28,40]. Methods which have an 
accurate energy and dipole moment (C and F) have good transport properties. 
Methods with a higher energy and a lower dipole moment (B and D) have trans- 
port properties which are too fast. The diffusion constant is larger and tnmr is 
smaller for these methods. For the model (E) with a lower energy and a larger 
dipole moment, the transport properties arc too slow. For method E, constant 
temperature dynamics is necessary to avoid heating. Transport properties are sen- 
sitive to how the velocity rescaling is done. It is preferable to use constant E,V,N 
dynamics but the Nose-Hoover method for constant temperature dynamics can re- 
produce diffusion properties well [41]. Constant temperature dynamics with the 
Nose-Hoover method were run using Ewald as a check and the resulting diffusion 
constant and tnmr were identical to the constant E,V,N results. In other studies 
with non-polarizable water potentials, the diffusion constant was found to be about 
the same as the Ewald result using spherical truncation [5] , a reaction field [5] and 
the shifting potential [4]. 

Structure. The radial distribution functions are sensitive to the treatment of 
the long-ranged interactions. Figure 2 shows the oxygen-oxygen radial distribution 
function, goo(r), for Ewald with 512 molecules and the shifted potential (method 
D) with 256 and 512 molecules. The Ewald goo(i')'s do not show a size dependence 
(data not shown) but the shifted potential results do show a size dependence. The 
first peaks do not have as much structure as Ewald, which is consistent with the 
smaller charges that result using the shifted potential. Also, there is structure 
around the cut-off distances (9.85 A and 12.4 A for the different sized systems) due 

TABLE 3. Dipole moment (Debye), translational diffusion constant 
(10~^ m^/s) and tnmr (ps). 

number of dipole moment diffusion constant tnmr 



molecules 


A 256 


2.62(1) 


2.0(3) 


2.2(2) 


512 


2.62(1) 


2.1(3) 


2.0(1) 


B 256 


2.51(1) 


3.1(3) 


1.1(2) 


512 


2.56(1) 


2.7(4) 


1.5(2) 


C 256 


2.62(1) 


1.9(2) 


2.1(1) 


D 256 


2.53(1) 


3.1(4) 


1.2(1) 


512 


2.55(1) 


2.7(3) 


1.4(2) 


E 256 


2.66(1) 


1.5(1) 


2.8(2) 


512 


2.65(1) 


1.8(1) 


2.6(1) 


F 256 


2.62(1) 


1.8(3) 


2.1(3) 


512 


2.62(1) 


1.9(1) 


2.2(1) 


experiment 




2.3" 


2.1'' 


a. Reference 

b. Reference 


38] 
39] 



to truncation effects. For non-polarizable models, the agreement between method 
D and Ewald is much better, although there still remains structure around the 

cut-off distance [3,4]. The results using method B are similar to the method D 
results. They show a similar size dependence, but do not have peaks at the cut-off 
distances. The interactions are sufficiently modified near the cut-off distance so 
that there are no truncation effects (see Figure 1). 

The results using the reaction field method agree well with the Ewald results, ex- 
cept for the peaks near the cut-off distance (Figure 3). The reaction field method 
also gives peaks at the cut-off distance for non-polarizable potentials [5]. With 
spherical truncation, the radial distribution functions are similar to the Ewald re- 
sults and are smooth at the cut-off distances (Figure 4). There is slightly more 
structure in the first peaks consistent with the larger charges. There is no system 
size dependence in the two spherical truncation results and the curves are indistin- 
guishable from each other. The goo(r) using method C is identical to the Ewald 
result and is not shown. 

From the radial distribution functions, the pressure results can be understood. 
The simulation methods which produce the largest pressure difference (models D 
and E) are also those which have peaks in the goo(r) at the cut-off region. These 
peaks will contribute to the virial and change the pressure. 

CONCLUSIONS 

For systems with polarizable charges, the long-ranged interactions are coupled to 
the charge distributions. Therefore, modifying the treatment of these interactions 
will modify the charges and the Coulomb interaction, q^qj/r, is changed not only by 
changing 1/r to S(r)/r, but also by changing q^ and q^. The results presented here 
indicate that the properties of liquid water are more sensitive to the treatment of 
the Coulomb interactions for polarizable systems than for non-polarizable systems. 
The reaction field correction to truncation (method F) is found to work well for the 
energy and transport properties, although it does give a large negative pressure. 
The difference between the reaction field method and shifting the potential by Sd{^) 
(Equation 1), which did not work nearly as well, are small (see Figure 1). Subtle 
differences in the treatment of the electrostatics can cause large differences in the 
results. Another truncation method which docs not work well, S(r)=erfc(Ar), can 
be made to give results almost identical to the Ewald results just by adding two 
constant terms which couple to the charges (compare methods B and C). These 
terms represent the self-term of Ewald (Equation 8) and a mean-field approximation 
to the Fourier-space term of Ewald (Equations 10 and 11). The fact that the mean- 
field approximation works so well suggests that fiuctuations in the Fourier-space 
and the forces from this term are not important. Once these terms are added back 
into model B, the charges become equal to the charges with the Ewald method and 
all the other properties including the energy, pressure, dynamics and structure are 
accurately reproduced. 
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FIGURE 1. Comparison of the truncation methods with the Coulomb interaction, 1/r (sohd 
line) showing erfc(Ar)/r (dashed line), the shifted potential S£)(r)/r (dot-dashed line) and the 
reaction field interaction (dotted line). 




FIGURE 2. Oxygen-oxygen radial distribution function with Ewald and 512 molecules (solid 
line), shifted potential, S£)(r), with 256 molecules (dashed line), and shifted potential with 512 
molecules (dotted line). 




FIGURE 3. Oxygen-oxygen radial distribution function with Ewald and 512 molecules (solid 
line), reaction field method with 256 molecules (dashed line), and reaction field method with 512 
molecules (dotted line). 




FIGURE 4. Oxygen-oxygen radial distribution function with Ewald and 512 molecules (solid 
line), spherical truncation with 256 molecules (dashed line), and spherical truncation with 512 
molecules (dotted line). 



